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Breakdown of the local density approximation in interacting systems of cold fermions 

in strongly anisotropic traps 
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We consider spin-polarized mixtures of cold fermionic atoms on the BEC side of the Feshbach 
resonance. We demonstrate that a strongly anisotropic confining potential can give rise to a double- 
peak structure in the axial distribution of the density difference and a polarization-dependent aspect 
ratio of the minority species. Both phenomena appear as a result of the breakdown of the local 
density approximation for the phase-separated regime. We speculate on the implications of our 
findings for the unitary regime. 



During the past few years, considerable experimental 
progress has been achieved in creating systems of strongly 
interacting ultracold fermions |l|. One of the primary 
motivations for this research effort is simulating strongly 
correlated electron systems such as high-temperature su- 
perconductors and quantum magnets. The number of 
particles used in experiments with cold atoms is large 
but not macroscopic, so they may also exhibit mesoscopic 
phenomena. In this paper, we discuss several impor- 
tant manifestations of such mesoscopic effects in spin- 
polarized mixtures of ultracold fermionic atoms in the 
vicinity of a Feshbach resonance [2, ISJ, U, la, l6| . Focus- 
ing on the BEC side of the resonance, we demonstrate 
that some of the unusual features observed in experi- 
ments by Partridge et al. [3|, including a double-peak 
structure in the axial distribution of the density differ- 
ence and a polarization-dependent aspect ratio of the mi- 
nority species [7| , appear as a result of strong anisotropy 
of a confining potential. Our starting point is a rigorous 
proof that neither of the two phenomena can take place 
when the confining potential is parabolic and the local 
density approximation (LDA) holds (this has also been 
noted in Ref. p]). A recent preprint by Zwierlein and 
Ketterle [8J argued that the unharmonicity of the confin- 
ing potential should have contributed to the double-peak 
structure observed in Ref. [3|. In this paper, we show 
that even for a harmonic potential, but for a strongly 
anisotropic trap, such as the one used in experiments of 
Partridge et al, beyond-LDA corrections give rise to both 
the double-peak structure in the axial density difference 
and a polarization-dependent aspect ratio. 

In this paper, we consider spin-polarized fermion mix- 
tures on the BEC side of the Feshbach resonance at 
T = 0. At low temperatures and deep into the BEC 
limit, all minority-component fermions are paired, form- 
ing stable bosonic molecules. So the system can be 
thought of as a Bose-Fermi mixture, where "bosons" 
are tightly bound molecules and "fermions" are excess 
or unpaired fermions of the majority species. Boson- 
boson (i.e., molecule-molecule) and boson- fermion (i.e., 
molecule-fermion) scattering lengths can be related to 



the scattering length between two species of fermions as 
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Before presenting a formal analysis, we provide a sim- 
ple qualitative picture of our results. The leading cor- 
rection to the LDA comes from including gradient terms 
in expressions for the chemical potentials of bosons and 
fermions. The gradient term for bosons corresponds to 
the kinetic energy term of the Gross-Pitaevskii equa- 
tion and has been thoroughly analyzed in the literature 
14| . Gradient terms for fermions have also been con- 
sidered in the context of nuclear physics (see Chap. 4 
of Ref. [15|). Interestingly, it turns out that gradient 
terms for fermions have a small numerical prefactor, and 
have a negligible effect in the considered region of the pa- 
rameters. For strongly anisotropic traps, gradient terms 
smooth the density distribution of bosons in the tightly 
confined radial direction at the boundary of the boson 
cloud (see the inset to Fig. [T]). This means that there 
are "extra" bosons in the radial direction compared to 
the LDA model. These "extra" bosons interact repul- 
sively with unpaired fermions and push them out in the 
axial direction toward the long tips of the trap, leading 
to a double-peak structure in the axial density of excess 
fermions (see Fig. [2|). We point out that we do not find 
a critical value of polarization beyond which the double 
peak structure appears, although the peak strength de- 
pends on polarization, scattering length, and the total 
number of particles in the system (see Fig. [3]). 

The evolution of the aspect ratio of the minority 
species is also easy to understand in the BEC limit that 
we consider. The aspect ratio of minority species is the 
aspect ratio of bosons (molecules). When the LDA ap- 
plies, the aspect ratio of bosons is equal to the ratio of 
confining frequencies (in experiments reported in Ref. [3| , 
this ratio is around 50). In the extreme limit of full polar- 
ization, one can consider a single boson (molecule) in an 
effective potential created by the confining potential and 
excess fermions. Solving a single-particle Schroedinger 
equation for the boson gives a wave function with the as- 



pect ratio equal to the square root of the ratio of confining 
frequencies (this corresponds to the extreme breakdown 
of the LDA, see also 12|). When the number of bosons 
is small but finite, one does not find such a dramatic de- 
crease of the aspect ratio, since even a small number of 
bosons leads to a change in the fermion distribution (see 
the inset to Fig. [2]). However, this argument explains 
a general trend of decreasing aspect ratio of minority 
species with increasing polarization. 

The microscopic model that we consider is the inter- 
acting Bose-Fermi Hamiltonian 
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Here $/ {^b) is a fermion (boson) operator and m is 
the original fermion mass. The interaction parameters 
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(121) with the parameters (H]) can be rigorously derived in 
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the dilute limit n^ a <C Ij't-/ a <C 1, where nb and 
Uf are local boson and fermion densities, from a single- 
channel model of a wide Fcshbach resonance [11] . It has 
been found out numerically 16] in the absence of density 
imbalance that fermion mixture is well described by the 
model of interacting bosons up to fc/a ^ 1. 

Within the LDA and assuming a harmonic confine- 
ment, the density profiles are given as solutions of the 
following equations: 

Hf (nb, Uf) + V{x,y,z) = consti, (3) 

libinb, Uf) + 2V{x, y,z) ^ const2, (4) 



where V{x,y, z) — 



2 / 2 , 2\ 



Notice that in 



the experiment of Ref . [3( , the trap is highly anisotropic 
A = lo^/lOz = 48.6 ^ 1. Boson and fermion densities 
depend on coordinates (a;, j/, z) through the potential 
V{x, y, z) only. Hence we should have identical densities 
of each species for points that have the same value of 
the confining potential. Thus the aspect ratio of both 
clouds should always be equal to the trap anisotropy A. 
Then 3D densities along the z axis, nb{z) = n(,(0, 0,z) 
and nf{z) = n^(0, 0,z), provide complete information 
about densities everywhere in the trap. From 3D densi- 
ties nbj{z) one can calculate axial densities as 



^(z) = dxdyuf vOr2+^2yy]Y2+P 
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For z > 0, the derivative of the axial density is 

dn1{z) ^ 2-K 



dz 



A2 



znf{z) < 0. 



(6) 



This provides a proof of the absence of peaks in the ax- 
ial density away from z = 0. We note that a similar 
statement can be made in the unitary and BCS regimes: 
if the density of the majority component is larger than 
that of the minority component at any point in the sys- 
tem, then the axial density cannot have a peak except at 
z — 0. Within the LDA and harmonic trapping, phase 
separation is signaled not by the appearance of a double- 
peak structure in the unpaired fermion axial density, but 
by the existence of an extended region where the axial 
density difference has a zero derivative. 

We now consider beyond-LDA corrections, which arise 
due to the spatial derivatives of the densities. At the 
mean- field level [l3|, correlations between the Bose and 
Fermi clouds are neglected and one treats the fermion 
(boson) density as providing an external potential for 
the bosons (fermions). Including the gradient corrections 
discussed earlier, up to second order in gradients of the 
densities, we find 



up to s 



E^ d 



■^gbb-nl 



2m V 



gbfUbUf 
107r2 



-Vir)inf 
36n/ 



2nb) 



In principle, for strong variations of the fermion density, 
higher-order terms in gradients have to be included. The 
gradient expansion above is not suitable for studies of 
shell structure |17| , but since these effects are small in the 
case we are studying, this level of approximation is suffi- 
cient. Even though p-wave superfluidity of fermions due 
to boson-induced interactions has been predicted (see, 
e.g., Ref. [5[), the value of the gap is exponentially small 
in 1/nJ a and cannot affect significantly the density pro- 
files obtained from the expression above. 

Taking variations of E with the density, one obtains 
corrections to the local chemical potential (valid in the 
region of nonzero fermion and boson densities). 
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Notice that beyond-LDA corrections to the fermionic 
chemical potential carry a small prefactor of 1/36 com- 
pared to the analogous term for bosons. 

To investigate numerically density profiles in the pres- 
ence of kinetic terms, we use the steepest descent method 
of functional minimization (used previously for bosons 
in anisotropic traps [18|). We investigate the effect of 
beyond-LDA corrections for typical parameters used in 
experiments (cf. Ref. [3]). In Fig. [l] we take harmonic 
confinement with ujz/2'k — 7.2 Hz and uj±/2t: = 350 Hz, 
a scattering length that corresponds to B = 754 G, num- 
ber of particles N = N^ + N-^ = 2Nb + Nf ^ 9.46 x 10*, 



n,(10 m") 




300 



400 



500 



600' 



FIG. 1: Rescaled 3D densities of fermions in the Bose- Fermi 
model ((2}. These correspond to density differences of the 
majority and minority components n/ = ni — ni. LDA 
result (dotted), nf{r — 0,z — x) (dashed), and nf{r — 
x/A, 2; = 0) (solid) for harmonic confinement ujz/'in = 7.2 Hz 
and ijj±/2n — 350 Hz, scattering length corresponding to 
B = 754 G, Ni + N^ = 9.46 x 10*, and P = 9.5%. Inset: 
Bose densities for the same conditions. 
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FIG. 2: Axial densities of fermions (n°j — n'^ — nj) within 
LDA (dotted) and in the presence of beyond-LDA corrections 
(solid) , for the same parameters as in Fig. [T] The inset shows 
the dependence of the anisotropy of the Bose cloud on polar- 
ization P, with other parameters being the same as in Fig. [T] 
Ra and Rr are obtained from the fits to the columnar density 
as discussed in the text. 



and polarization P = {N^ - Ni)/{Ni+N^) = Nfl{Nf + 
2Nb) ~ 10%. The LDA density profile is almost com- 
pletely phase separated due to the strong repulsion be- 
tween bosons and fermions that spatially overlap only in 
a small region. In Fig.[Tl we compare rescaled 3D density 
profiles in radial and axial directions with LDA results, 
and in Fig. [5] we do the same comparison for axial densi- 
ties and identical values of the system parameters. Mod- 
ification of 3D densities can be of the order of 30% for 
fermions, and is much stronger in the radial direction. 
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FIG. 3: Dependence of the strength of the double-peak struc- 
ture e on polarization P, with other parameters being the 
same as in Fig. [T] The inset shows how e changes with a when 
we move away from the point considered in Fig. [T] (marked 
with an arrow). Kf is defined as in Ref. §]. 



This is expected, since the gradient corrections in the 
radial direction are A^ « 2.4 x 10'^ times larger. We ob- 
serve a peak in the axial density, with the density at the 
maximum about 10% higher than at z = 0. The modifi- 
cation of the axial density is less prominent than changes 
in the 3D density profile due to the integration entering 
the definition of the axial density. Figure [3] shows the 
dependence of the relative strength of the double-peak 
structure (e) on polarization (P), where e is defined as 
e = [n'i(z)|maa; — n'^(0)]/n'^(0). It decreases with increas- 
ing P, since for a larger polarization a smaller fraction 
of the fermion cloud gets affected by the boundary of 
the boson cloud. Being essentially a finite-size effect, 
e also decreases with increasing total number of parti- 
cles. However, this dependence is pretty weak, since the 
size of the cloud depends weakly on the total number of 
atoms (only as N^'^ for the unitary regime within the 
LDA). The inset to Fig. [3] shows the dependence of e 
on scattering length (a). As a increases, the region of 
spatial overlap of the Bose and Fermi clouds shrinks due 
to stronger repulsion, and this leads to larger gradients 
and larger beyond-LDA corrections for the boson cloud. 
In addition, due to the enhanced repulsion, the fermions 
are more sensitive to the corrections to the boson cloud. 
Both of these effects lead to an increase of e for larger 
a. We present results as a function of Kfa, where Kf is 
defined as in Ref. |2,]. Notice, however, that there is no 
"universality" in this curve: indeed, increasing the total 
number for fixed a would increase Kfa, while e would 
decrease, as discussed earlier. Although Kfa in Fig. [3] is 
of the order of 1 , the main physics we are interested in 
takes place at the boundary of the boson cloud, where 
local kfa is considerably smaller, so the treatment of the 
system as a Bose-Fermi mixture is well justified. 

Beyond-LDA corrections increase the size of the boson 
cloud in the radial direction compared to the LDA result. 



thus reducing the aspect ratio of the bosonic cloud com- 
pared to A. As the polarization (P) increases, for a fixed 
total number of atoms, the size of the bosonic cloud de- 
creases and the importance of beyond-LDA effects grows. 
The inset to Fig. [2] shows the evolution of anisotropy as 
a function of polarization for harmonic confinement. Ra 
and Rr are defined as follows: (i) the columnar density 
7T.g(r, z) is obtained from nh{r,z), assuming radial sym- 
metry in the x-y plane; (ii) one then fits noninteracting- 
fermion density distributions to the columnar density 
along axial and radial directions, ng(r, 0) = A{1 — -^)^ 

and njj(0, z) — B{1 — ■^)^. Notice these functional forms 
are not always suitable and the fits provide just an es- 
timate of the radii. The aspect ratio of the fermionic 
cloud does not get significantly modified in a harmonic 
confinement. 

So far we discussed only the BEC regime, where con- 
trollable analytic theory is available. We developed a 
consistent qualitative picture of the appearance of peaks 
in the axial density difference, n'^ = n? — n?, and a de- 
crease of the aspect ratio of minority species with increas- 
ing polarization. Quantitatively, the numbers we obtain 
in the BEC regime are smaller than what is observed 
experimentally in Ref. [3| for the unitary regime, but 
of the same order of magnitude as the effects of unhar- 
monicity discussed in Ref. [8|. Now we will comment on 
the relevance of our findings for the unitary limit, where 
most experiments are performed. Results from the BEC 
regime should not be directly extrapolated toward the 
unitary regime, since the treatment of the system as a 
Bose-Fermi mixture starts to break down, and the form 
of beyond-LDA corrections in the unitary regime is not 
known. Hence the statements made in this paragraph 
are only speculative. One can estimate, within the LDA, 
the chemical potential in the center of the trap based 
on Fig. 2C of Ref. [3|: for a typical configuration it is 
~ 10?ia;_L, where ?iu!± is the radial level spacing. The 
chemical potential at the boundary of the inner cloud 
is about three to four times hio±, whereas the gradient 
contribution in the radial direction, which is neglected in 
the LDA, is ~ h'^A'^/[m{AR)'^] ~ hcj^ (where AR is the 
difference in axial radii). It is thus clear that beyond- 
LDA corrections might also be relevant for the unitary 
regime. The appearance of double peaks has been inter- 
preted in Ref. (3[ as evidence for phase separation of ex- 
cess fermions from a paired central core. Although peaks 
should not appear in the LDA, the picture we have in the 
BEC limit provides some support to this interpretation. 
When excess fermions spatially overlap with the paired 
core, they are less sensitive to beyond-LDA corrections, 
which are important only at the boundary of the paired 
cloud. As phase separation takes place, excess fermions 
are expelled to the boundary, where beyond-LDA effects 
become more important. These corrections are stronger 
in the radial direction, hence fermions are pushed in the 



axial direction to larger \z\. This may give rise to a pro- 
nounced double-peak structure in the axial density. The 
high anisotropy of the trap is important for this scenario. 

To summarize, we considered unbalanced fermions in 
the BEC hmit for T = 0. We proved a "theorem" 
that prohibits peaks in the axial density for z ^ 
within the LDA and for the harmonic trapping poten- 
tial. We showed that for strongly anisotropic confine- 
ment, beyond-LDA corrections produce a double-peak 
structure in the axial density, and change the aspect ra- 
tio of the inner cloud. We discussed the implications of 
our findings for the unitary regime. 
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